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ABSTRACT 

We simulate deep images from the Hubble Space Telescope (HST ) using semi-empirical models of 
galaxy formation with only a few basic assumptions and parameters. We project our simulations all 
the way to the observational domain, adding cosmological and instrumental effects to the images, and 
analyze them in the same way as real HST images (“forward modeling”). This is a powerful tool 
for testing and comparing galaxy evolution models, since it allows us to make unbiased comparisons 
between the predicted and observed distributions of galaxy properties, while automatically taking into 
account all relevant selection effects. 

Our semi-empirical models populate each dark matter halo with a galaxy of determined stellar mass 
and scale radius. We compute the luminosity and spectrum of each simulated galaxy from its evolving 
stellar mass using stellar population synthesis models. We calculate the intrinsic scatter in the stellar 
mass-halo mass relation that naturally results from enforcing a monotonically increasing stellar mass 
along the merger history of each halo. The simulated galaxy images are drawn from cutouts of real 
galaxies from the Sloan Digital Sky Survey, with sizes and fluxes rescaled to match those of the model 
galaxies. 

The distributions of galaxy luminosities, sizes, and surface brightnesses depend on the adjustable 
parameters in the models, and they agree well with observations for reasonable values of those pa¬ 
rameters. Measured galaxy magnitudes and sizes have significant magnitude-dependent biases, with 
both being underestimated near the magnitude detection limit. The fraction of galaxies detected and 
fraction of light detected also depend sensitively on the details of the model. 

Subject headings: galaxies: evolution -- galaxies: formation - galaxies: fundamental parameters — 
galaxies: general galaxies: statistics - large-scale structure of universe 


1. INTRODUCTION 

The Hubble Space Telescope (HST) has spent much 
of its operational lifetime staring into deep space, sur¬ 
veying galaxies in their infancy and youth. The moti¬ 
vation for these deep surveys (the Hubble Deep Field 
and its successors) was to obtain time-lapse images that 
would reveal how galaxies formed and evolved. While 
we have made great progress in interpreting the deep 
HST surveys, t his program remains challenging and far 
from complete dGalametz et, a, 1.1 [Midi : IGuo et all 1 201 .‘Hi 
iBouwens et al.ll20l4 and references therein). There are 
two major reasons for this. First, the samples of high- 
redshift galaxies have been severely edited by selection 
effects, primarily limits on flux and surface brightness, 
effectively biasing the observable universe toward bright, 
compact galaxies. Second, on the theoretical side, there 
are still many significant gaps in our understanding of 
the physical processes that affect the baryonic compo¬ 
nents of galaxies (stars, gas, and dust) and the radiation 
they emit. These uncertainties are reflected in the many 
free parameters of the semi-analytical models and in the 
analogous sub-grid physics of the hydrodynamical mod¬ 
els. In this paper, we present a new approach to the 
analysis and interpretation of deep galaxy surveys that 
addresses both of these issues. 

To account for selection effects, we create simulated 
HST images of model galaxy populations, and we then 
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analyze these images in the same way as real HST im¬ 
ages, extracting catalogs from the simulated images to 
detect and measure the fluxes, sizes, and other proper¬ 
ties of the galaxies. Our simulated images include both 
cosmological effects (projection along pencil beams, red- 
shifting of passbands, dimming of flux and surface bright¬ 
ness) and instrumental effects (point-spread function 
[PSF], pixelation, noise, sky backgrounds) for any given 
HST camera, filter, and exposure time. We then extract 
catalogs of objects in the images with the wid ely used 
SExtractor software (|Bertin fc Arnoutsl 1 19961) . Thus, 
our procedure automatically takes into account all rele¬ 
vant selection effects, allowing us to make unbiased com¬ 
parisons between the predicted and observed distribu¬ 
tions of luminosities, sizes, and other properties of galax¬ 
ies. This is the approach recommended by textbooks 
on statistical inference: map the predictions all the way 
into the observational domain and make the comparisons 
there, often called “forward modeling.” 

The forward modeling method in not sensitive to small 
errors in the luminosities, sizes, and other properties of 
galaxies, or even to the exact definitions of these quan¬ 
tities, because such errors affect measurements of both 
the simulated and real images in the same way. In other 
words, these errors cancel out of the comparisons of the 
simulated and real distributions of luminosity, size, and 
so forth. This is one of the main advantages of the for¬ 
ward modeling approach. 

In contrast, nearly all work in this field is based on 
the opposite, but simpler, approach of comparing predic¬ 
tions with observations in the theoretical domain (“back- 
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ward modeling”). Some exceptions are the reconstruc¬ 
tion of mock i mages or data starting from semi-analytical 


tion ot mock images or data starting from semi-an 
models, (e.g. iBlaizot et al.1 120051 : lOverzier et al 


or hydrodynamical simulatio ns le.g iLotz et al 


lDevriendtll2010t lMozenall2013H , although the latter suffer 
from unrealistic star formation histories and too rap id 
grow th of stellar masses at early times dBouche et al.1 
l2010h . In this paper, we go beyond the creation of 
mock galaxy images by deriving simulated distributions 
of galaxy properties and comparing them with the cor¬ 
responding observated distributions. 

To limit the number of assumptions and parame¬ 
ters in our models, we adopt a semi-empirical descrip¬ 
tion of galaxy evolution. This description is based on 
the evolution of dark matter halos in cosmological N- 
body simulations, which is now well understood, in con¬ 
trast to the evolution of the baryonic parts of galax¬ 
ies. The main assumption of the semi-empirical descrip¬ 
tion is that most of the information needed for simu¬ 
lating a population of galaxies is already encoded in 
the merger trees of their dark halos. Each halo is as¬ 
sumed to host one model galaxy, and the properties of 
that galaxy are then uniquely determined by those of 
its halo, including its mass and size. The advantage of 
this method is that it sidesteps much of the complex 
and uncertain baryonic processes in galaxy formation; 
the disadvantage is that it likely oversimplifies some as¬ 
pects of these processes. This semi-empirical description 
has been de veloped in numerous s t udies over the past 
decade (e.g.. lYale fc Ostriked l2004 jConrov et al.1 120061 


Conrov fe Wechslerl 120091 Guo et alt 20101: Moster et ahl 


20Hj Behroo^etaL 20101 


Guo fc Whitel 12014 


Moster et al.ll2013HKravtsovl 2015 : iReddick et al.l f2013h 


2013 


We derive the radiative spectrum of each simulated 
galaxy from its star formation history using stellar pop¬ 
ulation synthesis models. The star formation history in 
turn follows from the growth of its halo mass, including 
both smooth accretion and discrete mergers with other 
halos. The radiative spectrum also depends on the metal- 
licity of the stellar population and the absorption by gas 
and dust in the galaxy and by gas in the intergalactic 
medium. 

In our implementation of this method, we use cutout 
images of real galaxies in the Sloan Digital Sky Survey 
(SDSS) as templates for the visual appearance of our 
model galaxies, with their fluxes and sizes rescaled to re¬ 
flect galaxy evolution according to our models. This gen¬ 
erates much more realistic morphologies than the smooth 
Sersic bulge+disk light profiles commonly used in pre¬ 
vious semi-empirical or semi-analytical models and im¬ 
proves modeling of the detection incompleteness effects 
that arise from the internal dumpiness of real galaxies. 

Our approach gives us a valuable science tool for com¬ 
paring models of galaxy evolution. In this paper we build 
a reference model using plausible choices of parameters, 
and explore other models by changing one parameter at 
a time. Thus, we can test the sensitivity of the simu¬ 
lated universe to each of these parameters by compar¬ 
ing the statistics derived from their respective simulated 
images to those from reference model and real HST im¬ 
ages. Moreover, since the input properties of each model 
galaxy on the simulated images are known, we are able 
to quantify the model-dependent output galaxy detection 
efficiency by comparing the input and output distribu¬ 


tions of galaxy properties (such as luminosity or size). 
Our approach can be used to inform the design of fu¬ 
ture surveys (e.g., choice of filters and exposure times) 
by addressing directly the question of which data are 
most useful to discriminate between different theoretical 
models. This will be especially valuable in planning the 
deepest surveys with the James Webb Space Telescope 
( JWST ). 

We emphasize that the main purpose of this paper 
is to demonstrate the validity and utility of a general 
method for analyzing and interpreting deep galaxy sur¬ 
veys. This is an initial, exploratory study. We regard 
our specific implementation of the method and the first 
results obtained from it and presented here as being il¬ 
lustrative rather than definitive. There is scope for fur¬ 
ther development of the method and refinement of the 
results. Nevertheless, the overall agreement we find be¬ 
tween our simulations and observations— with no fine- 
tuning of parameters—is remarkable and encouraging. 

This paper is arranged as follows. fj2] explains the 
method for modeling and building simulated universes, 
from the selection of semi-empirical models and the dark 
matter simulation to the building of simulated HST im¬ 
ages. The main results, presented in (J3l include th e im - 
plementation of stellar mass-halo mass relations 1 33. ID . 
present-day mass and luminosity distributions derived 
from the semi-empirical models of a simulated universe 
(» and a comparison of the luminosity, size and sur¬ 
face brightness distributions extracted from the simu¬ 
lated images with measurements from real HST images 
( 33.3D . We devote §4 to an analysis of the cosmic star 
formation rate density. Lastly, §5] discusses the results 
and summarizes our conclusions. 

2. BUILDING A SIMULATED UNIVERSE 

In this section we describe in detail the steps required 
to build a self-consistent simulated universe and asso¬ 
ciated HST images. We also describe the parameters 
chosen for our reference model for the universe, as well 
as variations on those parameters explored in the other 
models. 


2.1. The Dark Matter Simulation 

Following standard practice, we use a ACDM simula¬ 
tion as the three-dimensional (3-D) skeleton of our sim¬ 
ulated universes, placing a model galaxy at the center 
of each dark (sub)halcQ. That defines the spatial distri¬ 
bution and number density of galaxies as a function of 
redshift. Most of the information needed in our method 
is in fact provided by the halo mass and size evolution 
along halo merger trees, which does not involve any fit¬ 
ting of free parameters (see 32.21 and 32.4( . 

For the merger trees, we use the milli- Millennium cos¬ 
molo g ical dark matter simulatio n (mMS) (Springc l et al.1 
120051 iLemson fc Springel 1200611 . Its smaller size com¬ 
pared to the full Mi l lennium or Millennium II simula- 
tions (iSoringel et al.l 120051 iBovlan-Kolchin et all 12009(1 
makes it easier to use for this exploratory work; we show 
below that the coarser mass resolution in the mMS is 
adequate for our simulations (' 33.3.111 . Halos are defined 
using friends-of-friends groups as explained in lGuo et al.1 

1 We u se “h alo” and “sub-halo” interchangeably, following 
I Guo et al.1 iFOlOH . 
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(|2010l) . Bound dark matter structures or (sub)halos are 
composed of the most massive main (or central) sub-halo 
surrounded by satellite sub-halos. 

The mMS has the same cosmology and particle mass 
(1.18 x 10 9 Mq) as the much larger Millennium simula¬ 
tion, but with both a smaller box size (85.62 Mpc) and 
a reduced number of particles (270 3 ). The cosmologi- 
cal parameters used are the ones obtained by WMAP 1 
(ISners-el et alJlihlilli . i.e., fl M = 0.2 5, Ha = 0 7 5. hn = 
0.73 and as = 0.9. As explained bv lGuo et al.1 (l2013aT) . 
the differ ence between the WM AP1 and the standard 
WMAP 7 (Komatsu et al. 1(201 ll ) cosmologies does not af¬ 
fect significantly the relevant aspects of dark matter 
structure. In fact, a smaller WMAP7 as = 0.807 is 
counterbalanced by a greater VLm = 0.272, which results, 
for example, in the WMAP1 and WMAP7-derived halo 
mass functions being very similar at z = 0. 

2.2. Constructing Stellar Mass-Halo Mass Relations 
with an Intrinsic Scatter 

We obtain the stellar mass of model galaxies in the sim¬ 
ulation using semi-empirical modeling. This approach 
defines the stellar mass M s of the galaxy as a function 
of the dark matter mass of the halo, M s = M s (Mh a io)> 
with the function M s (Mh a io) defined to be the stellar 
mass-halo mass (SMHM) relation. Although this is a 
simple one-to-one relation, it can be readily modified to 
includ e statisti cal scatter or redshift dependence (e.g., 
iBehroozi et all [20131 1. In this paper we adopt several 
SMHM relations from the literature, and use them for 
building our simulated images. We introduce a novel, 
self-consistent approach that naturally adds scatter to 
the SMHM mass relation. 

As a measure of the halo mass, we use the virial 
mass M v ; r (mass enclosed inside the maximum radius 
within which the mean density is 200 times the critical 
val ue), w hich i s ob tained from the value-added catalog 
of I Guo et al.1 (1201011 based on the milli-Millennium simu¬ 
lation. In this catalog, the mass of a central halo is given 
by M v i r , whereas for a satellite halo it is the maximum 
M v i r ever attained before becoming a satellite. In the 
senri-empirical approach, this preserves the stellar mass 
of satellite galaxies even while the outer parts of their 
dark halos are being tidally stripped as they orbit within 
a central halo. 

Figure [T| shows some of the SMHM relations found 
in the literature. Our reference model (or Model 1) 
for simulating the uni verse uses th e red shift-independent 
SMHM relation from I Guo et all ([2010f l. They computed 
this relation using a halo abundance matching technique 
that equates quantiles of the low reds hift stellar mass 
function from SDSS dLi fc Whitel [200 9) to those of the 
present-day halo virial mass function from the Millen¬ 
nium simulations. The quantile matching is made over a 
range of masses around the typical value M* (i.e., the 
“knee” in the SMHM relation as well as in the mass 
function), with extrapolations to the very low and high 
mass regimes, where both the stellar and halo mass dis¬ 
tributions are not well constrained. Anot her opt ion is 
the re dshift-evolving SMHM relation of Be hroozi et all 
(12013T) . This relation is based on stellar mass functions 
and cosmic star formation rates u p to z = 8. For o u r non - 
evolving Model 2, we adopt the IBehroozi et al.1 (|2013T ) 
SMHM relation evaluated at z = 0, while for our evolv- 


3 



Fig. 1.— Stellar mass M s versus halo mas s Mhaio re l ations 
test ed in this paper . Inclu ded are the models of I Guo et alii (2010) 
and IBehroozi et al.l (120131) , as well as a linear SMHM relation 
M s = 0.025M| li . l | o . The low and high mass tails are necessarily 
extrapolations of the fits due to the limited reso lution of the DM 
simulations (log Mhaio < 10. 8 for thelGuo et all (12010 ) model and 
log Mh a io < 10.3 for the IBehroozi et al.1 1201317 model). Note that 
this last model evolves with redshift, with M* (characterizing the 
knee in the SMHM relation) becoming smaller at higher redshift. 

ing Mode l 3 , we adopt the full redshift dependence of the 
IBehroozi et all SMHM relation. Our Model 4 does not 
involve halo abundance matching, but is a simple linear 
relation given by M s = 0.025Mh a i o , where the slope has 
been chosen by eye to coincide with the other relations 
in the vicinity of M*. We include this last model not 
because it is realistic but to study the sensitivity of our 
results to a SMHM relation that is very different from 
those of our firs t th ree models (based on the results of 
IGuo et akll2010l and IBehroozi et al.l[2Q13h 

In this paper, we do not explicitly impose random scat¬ 
ter in the SMHM relation. Instead we explore the scatter 
that emerges naturally from the dark matter simulation 
as a consequence of one basic assumption: we assume 
that the stellar mass along merger trees is a monotoni- 
cally increasing function of time. This is physically plau¬ 
sible because the stellar mass is concentrated in the cen¬ 
ter of halos due to dissipation in the baryons, and as a 
result it tends to be retained during mergers. This is the 
simplest physically motivated rule we have found for cre¬ 
ating consistent galaxy stellar masses from dark matter 
simulations. The alternative is to assume that galaxies 
lose and gain mass willy-nilly as halo masses decrease and 
increase and as sub-halos merge; that appears much less 
plausible based on our current understanding of galaxy 
formation and dynamics. Another alternative is to im¬ 
pose scatter on the SMHM relation in a predete rmined 
mann er (as in the approach adopted by Behrooz i et al.1 

(IMP 

Our assumption of monotonically increasing stellar 
mass naturally leads to scatter in the SMHM relation 
as a consequence of the following three related effects: 

1. In dark matter simulations, individual halos can 
decrease in mass from one time step to the next, 


























































4 


Taghizadeh-Popp et al. 


for example in events of tidal s tripping. I n that 
case, we follow the approach of 1 Guo et all (|2010T ) 
and do not reduce the stellar mass accordingly, but 
retain the stellar mass present before the decrease 
in M halo . 

2. In halo mergers, the dark matter mass of a de¬ 
scendant Mhaio,desc can be smaller than the sum 
of the progenitor dark matter masses SAfhaio.prog 
as some particles become unbound during the col¬ 
lision. This is not a rare occurrence in the simula¬ 
tions. In order to conserve the stellar mass content, 
we must break with the one-to-one fixed SMHM re¬ 
lation. 

3. Observed SMHM relations are intrinsically non¬ 
linear (see Figure [I]). That leads to a conflict 
with the assumed monotonic growth of M s in halo 
mergers. Consider the case when two halos with 
Mhaio,prog ~ M^o merge to create a new halo 
with Mhaio,desc > Mj( alo . The SMHM relation 
increases more slowly than linearly above M£ alo , 
which means that M s for the new halo is supposed 
to be smaller than the sum of the stellar masses 
for the merging halos. That conflicts with the as¬ 
sumption that the stellar mass along the merger 
tree must increase monotonically. 

To eliminate these conflicts we adjust stellar masses 
retroactively as follows. If a halo is found to have a stel¬ 
lar mass (according to the imposed SMHM relation) that 
decreases with time, we decrease the stellar masses of its 
immediate progenitors to enforce a monotonic increase in 
the stellar mass content of dark matter halos across their 
merger trees. This adjustment is applied recursively to 
all the progenitors of halos with modified stellar masses. 
Note that we are not assuming reverse causality with 
this scheme! Our assumption is that lower mass halos in 
denser environments (which are going to merge in the fu¬ 
ture) have their star formation rates suppressed by these 
environments. This procedure is described in more detail 
in Appendix lAl 

One natural consequence of this procedure is that, at 
a given halo mass, there is a scatter in the stellar masses 
that tends to fall slightly below the one-to-one SMHM 
relation at some points in the merger history. The results 
with the modified SMHM relations will be shown in EU 

2.3. Illuminating Galaxies in Dark Matter Halos 

The luminosity and spectrum of a model galaxy at any 
time are determined by its star formation history com¬ 
puted along the past merger history of its host dark mat¬ 
ter halo. Star formation is implied when the stellar mass 
of a descendant halo is greater than the total stellar mass 
of its progenitors in consecutive simulation time steps, or 
when a single halo increases its dark matter mass (and 
hence its stellar mass) due to accretion of surrounding 
dark matter particles. The stellar mass increase between 
simulation time steps implies a star formation rate, which 
is used in stellar population synthesis models to compute 
the emitted spectrum of the galaxy. We assume that the 
star formation rate between time steps is uniform, so that 
the star formation history is completely determined by 
the stellar mass history of a halo. Note that since we 



Fig. 2.— Probability density distributions of the (median sub¬ 
tracted) logarithms of R v \ r (dark matter halo virial radius) and R$o 
(SDSS galaxy r-band half-Petrosian-flux radius). Note that the 
radius distributions resemble Gaussian curves, as expected for ap¬ 
proximately log-normal random variables. The median values are 
(logf? v ; r [Mpc]) =-1.159 and (log /tgo[Mpc])=-2.824. with standard 
deviations oflogf? v ; r [Mpc]) = 0.176 and rrflog figofMpc]) = 0.285. 


have forced the stellar mass to increase monotonically 
with time, negative star formation rates are automati¬ 
cally excluded. 

We thus simply reconstruct the spectrum of the model 
galaxy (and derived photometry) as the sum of a series 
of uniform starbursts between each of the time steps in 
the simulation. When a galaxy is placed in a simulated 
image at a particular redshift z g (implying an age t g ), 
its rest-frame luminosity as a function of wavelength is 
computed using the star-formation history up to time 
t g using the evolutionary ste llar synthesis code GALAXEV 
bv lBruzual & Chariot] ( 2003 ). Note that the redshift z g 
is not required to be at one of the discrete simulation 
time steps, since we can integrate the star formation rate 
from the previous time step to the actual time t g . (The 
need of simil ar interpolation schemes has also been no¬ 
ticed by lYid 12010l ). In our reference model, we adopt 
the Chabrier initial mass function, with a fixed solar 
metallici ty (Z = 0.02) and the standard dust model from 
iCharlot fc Falll (120001) (■r„ — 1 and gL = 0.3). Our mod¬ 
eling of galaxy spectra is flexible, as we can in principle 
elaborate this model with additional variables, such as a 
redshift-dependent metallicity. 

2.4. Deriving Galaxy Sizes from. Dark Matter Halo Sizes 

We adjust the size of a model galaxy so that it is al¬ 
ways a fixed fr action of the evolving size of its dark mat - 
ter halo (e.g., iFall k. Efstathioul Il980t lKravtsovll2013D . 
We determine the constant of proportionality between 
the galaxy size and the halo size by comparing their re¬ 
spective z = 0 distributions, as in the halo abundance 
matching method. However, we match only their medi¬ 
ans instead of the full distributions, which suffer from 
incompleteness in the tails. It is well known that size 
distributions for halo s and galaxies are close to being 
log-normal (e.g., iShen et aLll2003l ). Since we choose to 
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make galaxy size linearly proportional to halo size, the 
scale parameter and the median subtracted logarithmic 
distributions of both distributions should be the same. 
Figure [2] shows the distributions of halo virial radii from 
Millennium and r-band half Petrosian flux radii R 50 for 
our main sample of SDSS galaxies from H2.51 The latter 
has been correcte d for incompleteness using the Vmax 
method (lSchmidtlll968i). as shown with a similar galaxy 
sample bv iTaghizadeh-Poop et aH (120121) . We find the 
relation i? 50 = 0.022I? v i r . The dispersions are different 
(indicating that our assumption is not completely accu¬ 
rate) but are similar enough for our purposes. As in the 
case of the spectra, we interpolate the galaxy size be¬ 
tween discrete time steps to the redshift z g of the model 
galaxy. 

2.5. Galaxy Image Cutouts from SDSS 

Our method places cutouts of real galaxy observations 
onto our simulated image, which has advantages over us¬ 
ing smooth analytic galaxy light profiles. Galaxies often 
have a clumpy internal structure, which affects source 
detections. A real galaxy could be detected as two or 
more separate objects, especially if its surface brightness 
approaches the image noise level. Real galaxy cutouts 
recreate this effect and are more realistic than smooth 
bulge and disk light profiles adopted in previous semi- 
empirical and semi-analytical models. Of course, this 
effect is not as important when the apparent galaxy sizes 
are comparable to the width of the point spread function, 
as may happen in the high redshift limit. 

A database of SDSS galaxy images is used for the 
cutouts. We select from the SDSS image database the 
real galaxy that is the closest neighbor to the model 
galaxy in a multi-dimensional space of galaxy proper¬ 
ties. Once the closest matching SDSS galaxy is found, 
we rescale its flux and size in order to make them equal to 
those of the model galaxy. No free parameters are fitted 
or required in this step. The full details of the SDSS data 
and the selection procedure is described in Appendix [Bl 

One might wonder whether SDSS galaxies are clumpy 
enough to be accurate models of high-redshift galaxies, 
since galaxies at higher redshifts tend to be more irreg¬ 
ular than local galaxies. The changes in morphology are 
due both to the shift of optical band filters into the rest- 
frame ultraviolet for distant objects and also to a higher 
merger rate and dynamically less-relaxed structures in 
the early universe. 

While these effects are worthy of further exploration 
in the future, for this paper the SDSS images are a good 
basis for simulations. Our analysis of galaxy counts and 
detection efficiencies relies on observations and simula¬ 
tions in the HST WFC3/IR F160W filter (A = 1.6 £tm); 
this implies that objects with redshifts less than 3 have 
an SDSS filter (from griz) that is at the appropriate rest- 
frame wavelength, and the different morphologies in the 
ultraviolet are irrelevant. Moreover, galaxies at redshifts 
beyond 3 tend to be sufficiently compact that their in¬ 
ternal dumpiness has little impact on their detectability. 
The median FWHM size of detected galaxies with z > 3 
in our simulations is 0.3 arcsec, which is only twice the 
FWHM of the 1.6 /im WFC3 PSF. While there is room 
for improvement in modeling the internal structure of 
galaxies, particularly for bluer filters in the 1 < z < 3 
redshift range, the SDSS cutout images are certainly an 


improvement over models that use smooth analytical pro¬ 
files. 

2.6. Assembling the Simulated Image 

Once model galaxy properties are calculated, we gen¬ 
erate pencil-beams through our simulated volume and 
project them onto the plane of the sky. We then simu¬ 
late ACS/WFC camera images, with their visible filters, 
as well as corresponding WFC3/IR images (sampled to 
the ACS/WFC pixel size), with their infrared filters. We 
include realistic PSFs for both cameras. 

To determine the 3-D structure within these pencil 
beams, we use a Monte Carlo method. This approach 
is much simpler and faster than other alternatives, such 
as the replication of a simulation box much smaller than 
the depth spanned by the simulated image. We first sam¬ 
ple a random redshifts z g from a distribution that gives 
constant comoving volume per redshift interval. Then we 
randomly select a dark matter halo at z = 0, choose all 
of its progenitors found at z g , and place them in the sim¬ 
ulated image (at redshift z g ) according to their relative 
3-D spatial positions in the simulation box (interpolat¬ 
ing between time steps). Although the large-scale cor¬ 
relations are discarded, we still preserve the short-range 
correlations between progenitors, while reducing consid¬ 
erably the computation time. 

The visual appearance of the model galaxy on the sim¬ 
ulated image is given by the best matching SDSS galaxy 
cutout (Appendix [B]). We select the SDSS band whose 
redshifted central wavelength A c (to redshift z g ) is the 
closest to the band of the simulated image. Due to the 
redshift of wavelengths, we end up using the bluest vis¬ 
ible SDSS bands to represent most high-redshift simu¬ 
lated galaxies in the visible and infrared HST filters, 
since SDSS lacks the UV and far UV bands that would 
be more suitable for this redshift regime@. Since the u 
band is noisy in SDSS, we use the g band as the bluest 
bandpass. In fact, all model galaxies placed on the simu¬ 
lated image are represented by a g-band cutout at z > 2 
for ACS/WFC images or z > 3 for WFC3/IR images. 

We rescale the flux of the image cutout to match that 
of the model galaxy. We apply to the model spectrum 
the effects of redshift, cosmological dist ance dimming and 
intergalactic absorption (as in iMadaul 119951) . before ap¬ 
plying the photometric filter response. 

We rescale the sizes of the galaxy cutouts placed on the 
simulated image according to one of the following rules: 

1. No size scaling: The proper size of the original 
SDSS physical galaxy is left intact, irrespective of 
the model galaxy size and hence the size of its 
halo. The physical size is scaled to the apparent 
size on the image using the angular diameter dis¬ 
tance Da{z 3 ). 

2. Size scaling: The proper size of the SDSS galaxy 
is scaled to match the physical size of the model 
galaxy, which in turn is a fixed fraction of the halo 
size (as described above). The apparent size is then 
computed from Da(z 9 ). 

3 A cross-match between SDSS and GALEX might be useful 
in the future for adding ultraviolet bands to our suite of galaxy 
cutouts, but that is out of the scope of this exploratory study and 
analysis. 
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TABLE 1 

Galaxy Evolution Models 


Model (short title) Details 

Model 1 • Stellar mass-halo mass relation f rom [ Guo et ah I (f20T0l) 

(Reference Model) • Dust model from [Uharlot Sz Falll d2000l) 

• Fixed solar metallicity (Z = 0.02) at all redshifts 

• Apparent sizes of SDSS galaxy cutouts on image are scaled to the theoretical size pre¬ 
dicted by the dark matter halo size 

• SDSS galaxy cutouts are chosen to be the closest match to the theoretically predicted 
model galaxy u — r color and stellar mass. 


Model 

2 

Same as Model 

Model 

3 

Same as Model 

Model 

4 

Same as Model 

Model 

5 

Same as Model 

Model 

6 

Same as Model 

Model 

7 

Same as Model 
for the angular 

Model 

8 

Same as Model 
SDSS to model 


1, but using the lBehroozi et al.l (120131) SMHM relation evaluated at z = 0. 
1, but using the lBehroozi et al.l (120131) redshift dependent SMHM relation. 
1, but using the a linear SMHM relation M s = 0.025Mj ia i o - 
1, but using no dust model for galaxies. 

1, but using a fixed very low metallicity (Z — 0.0001) at all redshift. 

1, but without rescaling the intrinsic size of SDSS galaxy cutouts (except 
diameter distance scaling, also applied in the reference model). 

1, but the PetroR50 radius as well as u — r and M s are used for matching 
galaxies. 


Note. — Description of eight different models used for building a simulated universe. Our reference model (Model 1) contains 
the most plausible parameters and sub-models. Other models are defined by changing one of these parameters at a time. 


In the first model, the z=0 size-mass relation of galax¬ 
ies is preserved at all redshifts, whereas in the sec¬ 
ond model, the size-mass relation evolves with redshift, 
driven by the growth of the dark halos. Galaxy sizes tend 
to be smaller at higher redshifts in the second model due 
to evolution in the halo sizes. 

To add instrumental effects, we convolve with the HST 
point spread function, apply the HST instrument ef¬ 
ficiency and pixelation, and add noise. The HST in¬ 
strument modeling uses the same software (pysynphot, 
Tiny Tim), instrumental parameters and sky background 
as used in the standard HST Exposure Time Calculators, 
providing a high fidelity model of the real HST perfor¬ 
mance. 

2.7. Source Detection and Photometry of Simulated 
Images 

The detection, extraction, and photometry of 
galaxies are per f ormed by running SExtractor 
ifBertin &; Arnoutsl 119961) on the simulated images. 
Here we dir ectly follow the metho d and parameters 
described in IGalametz et al.1 (1201311 and designed for 
the analysis of real HST images. The complete output 
catalog merges SExtractor runs using two different 
detection modes. The Cold mode is used for detecting 
bright and extended sources, while the Hot mode 
is optimized for extracting faint and small objects. 
After extraction, detected or output galaxies on the 
final simulated image are matched (using the detected 
position and luminosities) to the input galaxies on the 
same image before adding the instrumental effects. This 
provides a direct measure of the detection completeness 
by comparing what SExtractor detects to what was 
originally placed on the image. 

To compare our simulated galaxies to those from real 


galaxy surveys we use AUTO magnitudes and Petrosian 
radii i?p, as calculated by SExtractor. The Petrosian 
radii of our input galaxy cutouts, as given by the SDSS 
pipeline, differ from those returned by SExtractor on 
the already simulated images, probably due to different 
definitions for the radius in the SDSS pipeline and in 
SExtractor. Since we observe a linear offset between 
the distributions of these radii, we change the input 
SDSS values to match those from SExtractor by using 
logi?.p(SExtractor)=logi?p(SDSS) + 0.41. 


3. RESULTS 

The following subsections present our results for eight 
different models of simulated universes. These models 
are detailed in Table [T] where the reference model in¬ 
volves the most common choices of parameters as ex¬ 
plained in the previous sections, and the seven remaining 
models vary one parameter of the reference model at a 
time. In models 2-4 we change the SMHM relation; in 
models 5-8 we explore extreme values of other models 
parameter (often deliberately unrealistic values) to test 
the sensitivity of the observations to these parameters. 


























Simulating Hubble Images With Galaxy Formation Models 


7 




log(M ha i 0 [M 

sun]) 



log(M ha i 0 [M 

sun]) 



log(M halo [M 

sun]) 


log(M halo [M 

sun]) 


Fig. 3.— Modified stellar mass-halo mass relations for the four models tested in this paper, obtained after retroactively reducing the M s 
values predicted by the one-to-one SMHM relations in Fig. U to enforce a monotonically increasing M s as a function of time along merger 
trees. Data from three redshift time steps as well as the combination from all time steps is shown. The colors show the log-scaled number 
counts in bins of size AlogMhaio — 0.06 by A log M s = 0.01. 
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Fig. 4.— Cumulative distributions of the difference between the modified SMHM relation shown in Fig. [3] and the original one-to-one 
relation from Fig. \T\ (values from all redshift time steps are included). The 50%, 67% and 95% quantiles are indicated by dashed vertical 
lines. About half of the halos in each model have their masses unmodified (x = 0). 





































































































Fig. 5. — Present day stellar mass function for a simulated uni¬ 
verse according to the different models of stellar mass-ha lo mass 
rel ations tested in this paper: Model 1 (Guo ct ah 201(f), Model 
2 (IBehroozi et al.l[20130 and Model 4 (a linear SMHM relation). 
The open symbols show the stellar mass functions at z = 0 from 
the milli-Millennium and Millennium II simulations, computed by 
converting dark matter halo masses into stellar masses using the 
I Guo et al. (2010) relation. Obser ved stellar mass functions mea- 
sured from local SDSS galaxies dBaldrv et al.l [20081 ; ILi fe White! 
[2009h are also shown. 

3.1. Modified SMHM Relations with Natural Intrinsic 

Scatter 

As described in £12.21 and Appendix [A] we modified the 
fixed SMHM relation, selectively reducing the values of 
M s assigned to halos to enforce a monotonically increas¬ 
ing M s along merger trees. This leads to the natural 
dispersion of M s values shown in Figures [3] and Q] for 
the various SMHM relations explored in this paper. For 
all models, about half of halos do not have their stellar 
masses modified (meaning that they lie on the imposed 
SMHM relation). The scatter in the M s -Mh a i 0 distri¬ 
butions reaches about 0.07-0.12 dex in our simulations, 
similar to the scatt er inferred indirectly f rom theory and 
obser vations (e.g., iReddick et al.l 120131 IBehroozi et al.1 
120131) . The linear SMHM relation (Model 4) is the least 
affected by our modifications and shows the smallest dis¬ 
persion because it lacks the SMHM non-linearity present 
in the other models. For the non-linear Models 1-3, the 
scatter is smallest near Mh a io ~ 1O 12 M 0 because the stel¬ 
lar mass corrections are minimized near the peak of the 
SMHM relation. 

3.2. Simulated Universe Models at Low Redshift 

We evolve all models of simulated universes to the 
present day. In this section we compare the model prop¬ 
erties with low-redshift observations to test the overall 
accuracy of our method. 

3.2.1. Stellar Mass Functions 

Figure [3] compares the stellar mass function of model 
galaxies at 2 = 0 to those derived from SDSS. Th e ste llar 
mass function of the reference model, with its lGuo et al.1 


Fig. 6. — Global luminosity functions (r-band) calculated af¬ 
ter evolving our simulated universe to z = 0.1 under six different 
models. Also sho wn is the lumin osity f unction of local SDSS galax¬ 
ies measured bv IBlanton et ahl (|2005D . Note the similarity with 
Fig- El The downturn in the computed luminosity distributions at 
low luminosities is due to the finite particle mass resolution in the 
milli-Millennium simulation. 

(|2010f) SMHM relation, agrees w ith the observed stel¬ 
lar mass function from ILi fc White! ( 2009 ) in the range 
M s > 10 9 M q . That is expected since this SMHM re¬ 
lation was derived from halo abundance matching on 
mostly the same data. We also confirm that the stellar 
mass function predicted by the reference model fits per¬ 
fectly the stellar mass function derived by combining the 
Guo et al. SMHM relation with the halo mass function 
from the milli-Millennium simulation, as required by our 
semi-empirical modeling. However, the stellar mass func¬ 
tion of the reference model has an artificial downturn at 
M s < 10 'M 0 , which corresponds to the mass resolution 
of the milli-Millennium simulation in the identification 
of friends-of-friends groups (composed of a minimum of 
20 dark matter particles). Similar downturns are present 
in the stellar mass functions of the other SMHM rela- 
tions teste d. In contr ast, t he Millennium II simulation 
(|Bovlan-Kolchin et al.l f20Q9h . which has a dark matter 
particle mass 125 times smaller than the mMS, shows a 
power-law tail at the low-mass end. 

The halo abundance matching in the Guo model cov¬ 
ered the range 10 8 3 < M s / M 0 < 10 11,8 , which misses 
the upturn at M s < 1O 9 M 0 that is seen i n the more 
com plete SDSS stellar mass function from iBaldrv et al.1 
(2008) (sho wn by the black line in Fig. [5]) . The SMHM re¬ 
lation from IBehroozi et al.1 (120131) (Model 2, yellow line) 
fits this uptu rn much better since it was built using the 
IBaldrv et al.1 ( 2008 ) stellar mass function. 

The linear SMHM relation (Model 4) tracks the ap¬ 
proximately M~ 2 power-law mass distribution of dark 
matter halos. The proportionality constant for Model 4 
was chosen to roughly match the observed stellar mass 
function around its knee at M s ~ M*. Again, a clear 
mass resolution cutoff is present at lower stellar masses, 
also shown in the non-linear relations. This linear SMHM 
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model is obviously in strong conflict with the observa¬ 
tions at both higher and lower masses. 

3.2.2. Luminosity Functions 

Figure [6] compares the simulated and observed lumi¬ 
nosity functions at z = 0.1. To first order, the shapes are 
Schechter functions. The present-day luminosity func¬ 
tions from our models are similar in shape to their respec¬ 
tive stellar mass functions in Figure [5] This is expected 
since the r-band luminosity roughly traces the old stel¬ 
lar population that constitutes m ost o f the ste llar mass 
in present-day galaxies (e.g., iBell et ahl (1200311 1. Mod¬ 
els 5 (no dust) and 6 (low-metallicity) show good agree¬ 
ment with the iBlanton et al.1 (1200511 observed luminos¬ 
ity function in the high-luminosity tail. Since our refer¬ 
ence model contains dust and metals, its high-luminosity 
tail is shifted toward fainter luminosities with respect to 
those of Models 5 and 6. (Note that the luminosity in 
the r-band is reduced by dust absorption even though 
the bolometric luminosity is conserved.) On the other 
hand, the flatter slope at low luminosities for these three 
models mimics the flat tail present in their correspond¬ 
ing stellar mass functions, which is the result of applying 
halo ab undance matching to the stellar mass function of 
iLi fe White! (|2009H as explained earlier. A better fit is ob¬ 
tained f or Models 2 a nd 3 with the SMHM relation from 
iBehroozi et al.1 ( 2013 ). which includes the low-luminosity 
upturn. The luminosity function derived from the linear 
SMHM relation for model galaxies tracks the approxi¬ 
mately power-law mass distribution of dark matter ha¬ 
los. 

The perceptive reader may have noticed that even 
though semi-empirical modeling promises a perfect 
match of the predicted to observed universe, the global 
luminosity function predicted by our reference model 
does not agree perfectly with the observed luminosity 
function from SDSS at z ~ 0. T he reason fo r this is 
that the z ~ 0 SMHM relation from lGuo et al.l (12010(1 is 
not fully consistent with our method for converting stel¬ 
lar mass into light. They use the stellar mass function 
derived from ILi fo White! ( 200 9), which is not directly 
observed but is inferred from the observed SDSS galaxy 
luminosity function at z ~ 0 and an assumed mass-to- 
light ratio. On the other hand, our mass-to-light ratio 
is computed from th e dark matter halo masses and the 
SMHM relation using iBruzual fa Char loti (|2003f ) spectral 
models. We can in principle address this issue by com¬ 
puting our own SMHM relation using an iterative process 
that compares our measured luminosity functions to the 
observations. This level of precision is not needed in our 
present exploratory study, but it will be a useful longer- 
term goal for our approach. 

3.3. Simulated HST Images and Derived Statistics 

In this section, we show simulated HST images from 
our models and test the sensitivity of statistics derived 
from the images to the model parameters. We focus on 
the luminosity and size distributions, and we assess both 
biases in measured parameters and source detection in¬ 
completeness. For most tests we compare the perfectly 
known input values of size and luminosity (as given by 
the models) to the corresponding output values measured 
by SExtractor from the images. Our data comprise sim¬ 


ulated visible (ACS/WFC) and infrared (WFC3/IR) im¬ 
ages using filters and exposure times appropriate for the 
GOODS, CANDELS and HUDF surveys. Image sizes 
and pixel scales are those of a single ACS/WFC exposure 
(a 200" x 200" field with 0.0495" pixels). Note that this 
image area is small compared with the areas surveyed 
by many HST projects (e.g., GOODS and CANDELS), 
which encompass many ACS/WFC fields; we considered 
it appropriate to begin with a more modest sky area for 
this exploratory project. 

3.3.1. Results from the Reference Model 

Figure 0 shows a comparison between reference-model 
simulated and real HST ACS/WFC images. At first 
glance, the spatial distribution of galaxies and associ¬ 
ated sizes seem to be very similar. As expected, the 
HUDF-depth images show that many objects are hidden 
by noise in the GOODS-depth simulated images. 

Figure [8] shows scatter plots of the input (true model) 
values of the apparent sizes and surface brightnesses of 
both detected (using SExtractor) and undetected galax¬ 
ies in the HUDF-depth image. Most of the detected 
galaxies can be selected via a cut at tofi 60 W < 29, which 
does not depend strongly on redshift. In the lower panels 
of Figure [5J which plot surface brightness versus mag¬ 
nitude, it can be seen that there is a small, redshift- 
dependent population of low surface brightness galax¬ 
ies that are brighter than the magnitude threshold but 
nonetheless are not detected. 

Figure [9] plots the stellar masses of galaxies versus their 
input model apparent magnitudes. The stellar mass de¬ 
creases strongly with increasing redshift. However, an 
important conclusion from this plot is that even the least 
massive detected galaxies are still generally above the 
mass th resh old s set by both the SDSS survey (M s > 
10 6 M^. lKauffmann et al.ll2003f ) and the milli-Millennium 
simulation (M s > 10 71 M^. lGuo et aLll2010l) . Figure (TUI 
shows the typical apparent magnitude values found at 
the mass completeness limit of the milli-Millennium sim¬ 
ulation as a function of redshift. Although the median 
magnitude initially dims as redshift increases, the curve 
flattens out at mFieow ~ 32 for 2 > 3; although the 
galaxies at those redshifts have smaller masses, they also 
have higher star formation rates and younger popula¬ 
tions due to the requirement that their stellar masses be 
assembled in a short time. That approximately compen¬ 
sates for the cosmological dimming due to the increasing 
luminosity distance. 

Biases — As we noted above, Figure [8] shows that a 
relatively simple cut on the input (model) magnitude of 
galaxies predicts with reasonable accuracy which galaxies 
will be detectable in the images. However, detection is 
not the whole story. It is also necessary to measure the 
galaxy properties, and those measurements can be biased 
through complex effects related to the morphology and 
internal structure of the galaxy. Detection of an object 
does not ensure that its magnitude, colors, or size can be 
measured correctly. This is an area where our forward¬ 
modeling approach provides more reliable results than 
previous approaches. 

Figures |TT] and [12] show strong biases in the 
SExtractor-derived measurements of galaxy sizes and 
apparent magnitudes. From galaxies detected in the 
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Real GOODS image 



Simulated HUDF Image Real HUDF image 

Reference Model 


Fig. 7. — Simulated ACS/WFC F850LP+F606W+F435W images built from our reference model (GOODS and HUDF depths), compared 
to the equivalent real HST images. The same exposure times and display contrast are used for the comparison at each depth. The field of 
view is 2400 x 1200 ACS pixels 1/6 of the full ACS/WFC field of view). 


HUDF-depth simulated image, we compare the dif¬ 
ferences between the true (input) and SExtractor- 
measured ( output ) values of rriFi 60 W and logi?p. Our 
major conclusions are: 

1. There is a significant magnitude-dependent lumi¬ 
nosity bias, with the measured tofi6ow magnitudes 
on average fainter than the true galaxy magni¬ 
tudes (Fig.[lT]). For fainter galaxies, only the com¬ 
pact, high surface brightness nuclei are detectable 
above the image noise level, while the lower sur¬ 
face brightness extended envelopes are hidden in 
the noise. This bias is small for bright galaxies, 
but reaches median values of ~ 0.4 magnitude and 
extreme values of > 1 magnitude at the detection 
limit. This effect does not depend strongly on red- 
shift. 

2. There is also a strong bias in the measured sizes 
near the magnitude detection limit (Fig. fill) . 
The SExtr actor-measured sizes are smaller than 
the original input sizes of model galaxies around 
mFi 60 W — 29. The explanation appears to be the 
same as for the luminosity bias: the detectable part 
of a faint galaxy is significantly more compact than 
the true Petrosian radius due to masking of the 
more extended component by noise. The differ¬ 
ence between output and input values of log Rp 
has a median of 0.15 dex, with maximum values 
reaching even to 0.8 dex. Note that there is also a 
Malmquist bias affecting these distributions: faint 


extended objects may be detectable only if noise 
fluctuations make their surface brightness appear 
brighter than the detection limit. 

Figure |TT| also reveals a small bias toward larger mea¬ 
sured output sizes for sources ~ 1 mag brighter than the 
detection limit. This bias is largest in the z > 6 red- 
shift bin, where the difference reaches around 0.2 dex, 
and might be related to intrinsic peculiarities in the way 
SExtractor measures sizes. 

These effects are certainly detection biases and not 
a trend resulting from supposed evolutionary effects in 
the models. Figure [12] shows a scatter plot of the mea¬ 
sured (output) sizes versus measured magnitude for both 
GOODS and HUDF-depth simulated images. The strong 
bias in the SExtractor-measured sizes at the magni¬ 
tude detection limit are easily visible at both depths, 
and the bias begins at a magnitude that is determined 
by the detection limit rather than at any physically de¬ 
termined luminosity. Note that the underlying images 
(before adding noise) are exactly the same in these two 
cases: they are equivalent to observing the same sky re¬ 
gion twice with different exposure times. This rules out 
any artifact introduced by the model. Note how the mea¬ 
sured Rp extends down to values close to the PSF size 
(FWHM=0.151" for the F160W band image) at the mag¬ 
nitude detection limit. 

We note that the true biases could be even larger than 
those measured in our simulations. Real galaxy light 
profiles have extended wings, whereas our galaxy image 
cutouts only include the light present within two Pet- 
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Fig. 8. — Top : Log Petrosian radius versus apparent F160W 
magnitude for the HUDF-depth simulated HST image derived from 
our reference model (Model 1), separated in four redshift bins. 
Blue points show galaxies detected by SExtractor; orange points 
show undetected galaxies. For all galaxies, the input (model) radii 
and magnitudes are shown rather than values measured from the 
images. The black line is the same in all four panels and shows 
the median of the log Rp versus magnitude distribution integrated 
over all redshifts. Bottom : Apparent surface brightness p versus 
magnitude using the same color coding. 

rosian radii. According to iStrauss et al.l (l2002h . 82% 
of the light from a de Vaucouleurs profile is contained 
within this radius, and 99% of light from an exponen¬ 
tial profile is included. In addition to these moderate 
effects, there might be cases when the outer wings of 
the light profile are much more extended, with surface 
brightness falling close to R~ 2 . Then the light in the 
extended halo would dominate the total luminosity, but 
such extended emission would certainly not be detectable 
for faint galaxies, and the biases would be increased. The 
true radial dependence of the outer wings of galaxies is 
not currently well constrained, so the amount of light 
missing from these wings remains an open issue. 
Detection Efficiencies — The detection efficiency is an 
important statistic, since it tells us the amount of under¬ 
lying information that we are losing when observing with 



Fig. 9. — Stellar mass M s versus apparent F160W magnitude for 
the HUDF-depth simulated HST image derived from our reference 
model (Model 1), separated in four redshift bins. Color coding is 
the same as in Fig. [8] The lower dashed line s hows the smallest 
stellar mass for SDSS galaxies in the catalog of IKauffmann et al.l 
(|2003D (M s = 10 6 Mq). The upper dashed line shows the stellar 
mass threshold for halos resolv ed by the m ill i-M illennium simula¬ 
tion (M s = 10 71 Mq) using the I Guo et al.l (120101) SMHM relation. 
Note that images from our reference model observed at the HUDF 
depth are not affected by either threshold, since galaxies located 
below the thresholds are too faint to detect. 



Fig. 10.— Median apparent F160W magnitude as a function of 
redshift using our reference model for galaxies with stellar masses 
near the milli-Millennium stellar mass limit of M s = 10 71 Mq (up¬ 
per dashed line in Fig. EJ. The curve flattens at z > 3 where 
the youth and high star-formation rates of galaxies compensate for 
cosmological dimming. 


HST. In this paper, we consider two different measures 
of the efficiency: the number count detection efficiency is 
the fraction of galaxies on the image that are detected in 
the simulated data, and the light detection efficiency is 
the fraction of the total galaxy flux in the observing band 
that is detected. These efficiencies are easily computed 
from the model images since we know the properties of 
both the detected and undetected galaxies. Both quan- 
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Fig. 11.— Top : Magnitude biases: the difference between the 
true (input) mpigow magnitude of galaxies and the SExtractor- 
measured output magnitude as a function of the input magnitude. 
The results for four different redshift ranges in HUDF-depth sim¬ 
ulations are shown. Blue lines show the median magnitude dif¬ 
ferences over all redshifts and are the same in each panel. Points 
below the dashed line have measured magnitudes that are fainter 
than the true magnitudes. Fluxes are systematically underesti¬ 
mated for fainter galaxies. Bottom-. Size biases: same as top, 
but showing the difference between the true input log Up sizes 
and the SExtractor-measured sizes. There is a slight bias toward 
larger sizes for brighter galaxies, while fainter galaxies have mea¬ 
sured sizes that are substantially underestimated. Only the high- 
luminosity cores of faint galaxies are detected, while the remaining 
extended emission is lost beneath the noise. 


tities can be computed as a function of other parameters 
such as redshift or apparent magnitude. The number 
count efficiency is more commonly used but can be ill- 
defined when one considers the possibility of numerous 
faint galaxies that are undetected but contribute little 
stellar mass or light. The light detection efficiency is bet¬ 
ter behaved, since the integral of light from many faint 
galaxies is typically a small correction to the light from 
objects near the knee of the galaxy luminosity function. 

The top panel of Figure [T5] shows the detection effi¬ 
ciencies for number counts and for F160W-band light at 
HUDF depth as a function of the model input magnitude 
of galaxies. Both efficiencies drop sharply to zero at the 
magnitude detection limit of TOfi 60 W — 29. Interest- 




Fig. 12.— Measured (output) size log Up as a function of mea¬ 
sured magnitude mFl 60 W for detected galaxies. The left and right 
panels show results from simulated images at GOODS and HUDF 
depth, respectively. The blue lines show the median trend. Note 
that the tendency to underestimate galaxy sizes at the magnitude 
completeness limit is evident at both depths, which demonstrates 
that it is not the product of any peculiarity in the reference model. 


ingly, the number count efficiency has the same shape in 
all redshift bins, showing a slight decline toward fainter 
magnitudes but remaining above 80% before reaching the 
detection limit. It reaches values closer to unity only at 
bright magnitudes in the z < 2 redshift bin. The light de¬ 
tection efficiency has a similar behavior but lower values: 
it falls to ~ 60%-70% just above the magnitude com¬ 
pleteness limit. The lower efficiency for detecting F160W 
light is due to the magnitude bias discussed above: not 
only are galaxies undercounted near the detection limit, 
but the fluxes of detected galaxies near the detection 
limit are also significantly underestimated (Fig. fTIl). 

The bottom panel of Figure [T51 shows the detection ef¬ 
ficiencies as a function of redshift. These curves are a bit 
more complicated to interpret because they are affected 
by the finite mass resolution of the milli-Millennium sim¬ 
ulation. Very low-mass galaxies (< 10 '-'Mq) cannot ap¬ 
pear in the simulations, so their contributions to the un¬ 
detected counts and luminosity are omitted. That means 
that the computed efficiencies are only upper limits to 
the real values. This has little effect on the efficiency as 
a function of magnitude since very few low-mass objects 
would be brighter than 29th magnitude (Figs. [9] and [10]), 
but it has a significant effect on the redshift-dependent 
efficiencies since faint galaxies can appear at all redshifts. 

Figure fldl shows two attempts to correct the detection 
efficiencies for the low-mass galaxies. The solid curve 
simply omits all galaxies fainter than the median mag¬ 
nitude of the lowest mass halos in the milli-Millennium 
simulation (Fig. ITOl) : the efficiencies plotted are then ef¬ 
fectively the fraction of galaxies or light from galaxies 
brighter than ~ 32 mag, which is an upper limit to the 
true detection efficiency. The dashed curve instead inte¬ 
grates the light in the extrapolated power-law tail of the 
faint galaxy luminosity function (Fig. 1141) and includes 
the light of the missing faint galaxies as part of the total 
flux. This produces a modest correction in the light de¬ 
tection efficiency. Note that a version with extrapolated 
number counts is not shown because the number count 
correction factor is much larger (and far more uncertain). 

The number count efficiency decreases rapidly from 
80% at low redshift to 20% at z ~ 7. On the other 
hand, the light detection efficiency drops more slowly, 
reaching values of ~ 50%-70% at z ~ 7. The light detec¬ 
tion efficiency is greater than the number count efficiency 
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Fig. 13.— Top: SExtractor detection efficiency as a function 
of input raFi 60 W? measured on a simulated HST image at HUDF 
depth and built from our reference model (Model 1). The num¬ 
ber count detection efficiency is the ratio of the number of galaxies 
detected by SExtractor to that of all the galaxies placed in the im¬ 
age. The light detection efficiency is the ratio of the output fluxes 
(F160W band) of detected galaxies as measured by SExtractor to 
the corresponding model input fluxes of all galaxies placed in the 
image. Note that light efficiency is smaller than that for number 
counts because galaxies are not only missed but also have their 
fluxes underestimated. Bottom: Efficiencies as a function of red- 
shift. The solid lines include only galaxies brighter than the me¬ 
dian mpi 60 W ^ the stellar mass incompleteness limit (Fig. [lOj. 
The dashed line uses a power-law extrapolation bey ond the mass 
limit to estimate the light of faint galaxies (Fig. [TDl . 


because every redshift bin includes numerous faint galax- 
ies, and the detected ones (at fixed redshift) are generally 
the brightest and carry most of the light content of the 
image. 

The correction from including the extrapolated light 
from faint low-mass galaxies increases at higher redshifts. 
This is can be seen directly in the steepening of the slope 
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Fig. 14.— Apparent magnitude distribution of our reference 
model for HUDF depth separated in redshift bins. The green line 
shows a power-law low-luminosity tail fitted with slope a. The 
power-law slope at the faint end tail becomes steeper with increas¬ 
ing redshift. 


of the fitted power-law tail at increasing redshifts which 
changes from a ~ —1 at z < 1 to a = —1.75 at z = 7 — 8 
(Fig.[lT]). Interestingly, the steep slopes for young galax¬ 
ies at high redsh ifts are also found in the present-day 
universe. iTaghizadeh-Popp et all ( 2012 ) isolated similar 
populations of small, blue galaxies with rapid star forma¬ 
tion at z = 0 and also found faint-end luminosity slopes 
close to a ~ —1.6. 

Luminosity and size distributions — The luminosity and 
size distributions for all galaxies detected on simulated 
images from the reference model are shown in Figure fl5l 
together with the respective distribution of model input 
values and SExtractor-measured output values. When 
compared to the observed luminosities and sizes derived 
from the CANDELS GOODS-S Multi-wavelength Cat¬ 
alog (IGuo et al.ll2013bl) . the SExtractor output values 
from the reference model (red) agree surprisingly well. 
For the luminosity distribution, a good fit is found at in¬ 
termediate magnitudes, but the model predicts slightly 
too few galaxies at the bright end (possibly affected 
by sample variance due to the small simulated image 
area). The drop in the input luminosity distribution at 
WF 160 W ~ 32 is an artifact of the finite mass resolu¬ 
tion in the dark matter simulation. Since the detection 
limit of the HUDF-depth image (tofi 60 W ~ 30) is 2 mag 
brighter, we are safe from this artificial incompleteness. 
The SExtractor size distribution also agrees with the 
observations, although the model seems to be slightly 
shifted toward smaller sizes. The bias between the input 
and output distributions of detected galaxies is also vis¬ 
ible in this figure. The HUDF-depth image shows that 
the peak value of the output I?p distribution is shifted 
toward bigger radii by ~ 0.1 dex with respect to the 
input distribution. Since this bias is not as strong in 
the GOODS-depth image, we conclude that the ampli- 
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Fig. 15. — Galaxy number counts as a function of apparent 
F160W magnitude (top) and Petrosian radius (bottom) from sim¬ 
ulated HST images at both GOODS (left) and HUDF (right) 
depth using our reference model. The input values of all model 
galaxies placed in the simulated image are shown as orange 
squares. The input values for galaxies detected by SExtractor 
are shown as blue triangles, while the output values as measured 
by SExtractor for those detected galaxies are shown as red cir¬ 
cles. SExtractor-derived distributions are not corrected for in¬ 
completeness. Black lines show the actual observed distributions 
calculated from the CANDELS GOODS-S Multi-wavelength Cat¬ 
alog dGuo et al.H2013bh . There is good agreement between the 
simulated and observed distributions except for the smallest size 
objects. 


tude of the bias depends on the properties of galaxies 
near the detection limit, which differs in the GOODS 
and HUDF-depth images. Another interesting feature is 
that the output size distribution departs from the input 
size distribution at small values of Rp (down to the PSF 
size), where measured sizes for small galaxies are much 
less than their true sizes. This effect was also discussed 
above CFig. fl2l). 

3.3.2. Results from Other Models 

The statistics derived from simulated images built with 
other models are qualitatively similar to those of the 
reference model and follow similar trends. However, 
the simulated universe is sensitive to changes in the 
model parameters, and we can easily distinguish vari¬ 
ations among the derived statistics for different models. 

Figure [16] shows the simulated images from all models. 
At first glance it is easy to recognize differences between 
the reference model and the others. For example, some 
images show a deficit of small galaxies (e.g., Models 2 
and 3) or an excess of them (e.g., Model 4). Models 5 
and 6 with no dust and low metallicity look brighter, and 
Model 7 with no size scaling has much larger galaxies. 

Figure |T7] shows the size versus apparent magnitude 
for two models that stand out from the others. Model 
4, with the linear SMHM relation, predicts higher stel¬ 
lar masses in small-mass halos than the reference Model 


1 , especially in the high-redshift universe where halos 
are just starting to assemble. The net effect is to make 
compact galaxies brighter. Compared with the reference 
model in Figure [5] a much higher fraction of the galax¬ 
ies are detected by SExtractor, especially in the z < 2 
redshift range, since most of the galaxies on the image 
are brighter than the magnitude detection limit. Model 
7, with a non-evolving size-mass relation, has a very 
different distribution of apparent sizes compared with 
that of the reference Model 1 and the observed distri¬ 
bution. This is a powerful demonstration that galaxies 
were smaller in the past, even at fixed mass. 

The median apparent magnitudes at the mass com¬ 
pleteness limit have different redshift dependences among 
the models. In Figure [TS] most of the models follow the 
same behavior as Model 1, with the exception of Models 4 
and 5. Model 4 has median magnitudes at the mass com¬ 
pleteness limit 1.5-3 mag brighter than those of Model 
1, since the linear SMHM relation assigns more stellar 
mass content per unit halo dark matter mass. Model 
5, with no dust, brightens at 4 < z < 7 in comparison 
with Model 1, indicating that dust typically dims galax¬ 
ies about by ~ 1 mag in the redshifted filter bandpass. 
All the models, however, flatten at high redshift as dis¬ 
cussed for the reference model. 

Figure [12] shows detection efficiencies for all models 
as a function of both magnitude and size. The number 
count efficiency as a function of magnitude varies for all 
models similarly to that of the reference model, declining 
slowly to 80% just above the magnitude detection limit 
and then dropping dramatically, irrespective of redshift. 
The light detection efficiency as a function of magnitude 
also shows a close similarity among the models, dropping 
to 60%-70% just above the magnitude detection limit. 

Regarding the number count detection efficiency as a 
function of size, most models have high efficiency (80%- 
90%) for larger objects with a sharp drop at input radii 
of ~ 1". The drop is due to smaller objects also hav¬ 
ing smaller fluxes so that they are closer to the detection 
limit. There are, however, two models that stand apart. 
The linear SMHM relation in Model 4 makes compact, 
low-mass galaxies that are more massive and brighter 
than those in Model 1, and this greatly enhances the 
detectability of smaller objects. As a consequence, the 
detection efficiency stays above 50% at input sizes down 
even to 0.1" for z > 4. At the other extreme, Model 7 
(with the non-evolving size-mass distribution) has larger 
sizes and much lower detection efficiency at high redshift 
compared to the other models. The larger size distribu¬ 
tion is seen because galaxy sizes are not scaled in pro¬ 
portional to the halo sizes. Consequently high redshift 
galaxies are large fuzzy objects that are strongly sup¬ 
pressed in the images by the (1 +z) A cosmological surface 
brightness dimming, and very few of them (10%-20%) 
are detected. The smaller sizes of galaxies at high red¬ 
shift in most of the models are essential to making them 
detectable. The light detection efficiency as a function 
of size shows similar effects, although the differences for 
Models 4 and 7 are less dramatic. 

The detection efficiencies as a function of redshift in 
Figure 1^01 are only upper limits to the real values, since 
we are missing halos smaller than the mass completeness 
limit in the dark matter simulation. We calculate these 
efficiencies using only galaxies that are above the median 
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Fig. 16.— Simulated ACS/WFC F850LP+F606W+F435W images (HUDF depth) derived from for all models explored in this paper. 
The same exposure times and display contrast are used for the comparison at each depth. The field of view is 1200 X 1200 ACS pixels 
(~ 1/12 of the full ACS/WFC field of view). Changes in the distribution of luminosity and sizes are apparent for different models. 


magnitude set by the simulation mass resolution. The 
distributions for different models are qualitatively simi¬ 
lar although there are some significant variations. The 
efficiency amplitudes depend strongly on the model and 
on how compact and luminous the galaxies are, especially 
in the case of the number count efficiency. For example, 
Model 4, with its more compact galaxies, shows by far 
the highest number count detection efficiencies, > 60% at 
all redshifts. The light detection efficiency is higher than 
the number count efficiency for all models. It is notewor¬ 
thy that Model 5, with no dust, shows the highest light 
detection efficiencies at z > 4. These experiments show 
that the efficiencies depend on details of the models and 
that understanding them is key to inferring the physics 
of galaxy formation from observations. 

Figure [Til shows luminosity functions, which are espe¬ 
cially suitable for exploring the effects of our different 
models on the simulated images. Changing the SMHM 
relation modifies the stellar mass distribution across red- 
shift, which can be seen directly as a change in the appar¬ 
ent luminosity functions. According to Figure 0 the lin¬ 
ear SMHM relation (Model 4) predicts a much higher M s 
at a given Mh a i 0 than the Guo et al. relation, especially 
for halos with Mhaio /5 10 11 Mq. That shifts the peak 
of the luminosity function toward brighter magnitudes 
compared with the reference model. Since the galaxies 
are brighter, more of them are detected by SExtractor 
or seen by eye (e.g., Figs. HblorUTl). With the Behroozi 
et al. SMHM relation, M s is higher than with the Guo et 
al. SMHM relation at small values of Mhaio, but lower at 
high Mhaio (for any redshift). The net effect is a shift of 
the peak to brighter magnitudes, but at the expense of 
fewer galaxies just brighter than the peak. At z = 0, the 


Behroozi et al. SMHM relation is similar to that of Guo 
et al., so the luminosity function of Model 2 is closer in 
shape to that of Model 1. Galaxies in Model 5 (with no 
dust) are, as expected, more luminous than in the refer¬ 
ence model, in some cases ~ 0.7 magnitudes brighter in 
the range 22 < mFi 60 w < 25. Model 6 (with low metal- 
licity) behaves similarly to Model 5, but the magnitude 
shift is not as large. These 2 models also give good fits to 
the observations. The last 2 models share the same input 
z = 0 luminosity function as the reference model, since 
only the size of galaxies are modified. Therefore, Model 
8, which uses R 50 in the matching between model and 
SDSS galaxies, has a luminosity function at z = 0 similar 
to that of Model 1. However, Model 7 is slightly different 
at the onset of incompleteness, since here galaxy sizes are 
not scaled and are therefore larger than in the reference 
model. The galaxy sizes affect the surface brightnesses 
and so are influential in determining the detectability of 
faint galaxies. 

The size distributions in Figure [22] mostly behave in 
the same way as in the reference model, except for Mod¬ 
els 4 and 7. For Model 4 (linear SMHM relation), the 
peak of the output size distribution shows a larger bias 
toward bigger sizes (by ~ 0.2 dex) compared to the input 
distribution. Note also that the output size distribution 
for smaller detected galaxies always falls below that of 
the input distribution. The input sizes of the smallest 
detected galaxies extend well below the FWHM of the 
PSF (0.151"), but those bright enough to be detected 
have SExtractor sizes comparable to the FWHM. (We 
have not modified the SExtractor sizes to remove the 
effects of the PSF.) The size distribution for Model 7 is, 
not surprisingly, very different from the other distribu- 
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Fig. 17. — Petrosian radius versus apparent F160W magnitude 
for the HUDF-depth simulated HST image derived from Models 
4 (top) and 7 (bottom), separated in four redshift bins. Colors 
are the same as Fig. [8] The solid line shows the median values 
for these models, while the dashed line shows the median for the 
reference model for comparison. Small galaxies are given a higher 
luminosity in Model 4; the large galaxy sizes in Model 7 lead to 
many galaxies being undetected due to low surface brightnesses. 

tions, since its galaxies lie on a non-evolving size-mass 
relation rather than having their sizes rescaled by the 
evolving sizes of their halos. That results in many large, 
faint galaxies at high redshift that are difficult to detect. 

4. COSMIC STAR FORMATION RATES 

The cosmic star formation rate density (SFRD, the 
mass of stars formed per unit time and per unit co¬ 
moving volume) and its evolution with redshift have 
played a key r ole in studies of galaxy formation (see 
iMadau fe DickmsonI[20141 and references therein). We 
have postponed consideration of the SFRD until now 
because, although this relation is derived from observa¬ 
tions and is often portrayed as an “observed” relation, 
it is unsuitable for direct com parison with_our models. 
In previous studies (e.g., iBehroozi et al.1120131) . predic¬ 
tions of the SFRD are compared with observations in the 
theoretical domain, i.e., by combining observed galaxy 



Fig. 18.— Median apparent F160W magnitude of galaxies hav¬ 
ing stellar masses near the simulation mass incompleteness limit 
of logM a /MQ = 7.1, shown for galaxies in different redshift bins. 
Data is shown for all the models explored in this paper. 

counts and colors with corrections for dust absorption 
but without corrections for missing light at low lumi¬ 
nosity and surface brightness. We argue instead that it 
is better to project theoretical models into the observa¬ 
tional domain by creating simulated data and making 
comparisons directly with the observed quantities (“for¬ 
ward modeling”). However, a complete study using that 
approach is beyond the scope of this paper. 

Nonetheless, a comparison between our model proper¬ 
ties and SFRD values from the literature is illuminating. 
Figure l2dl' a) shows the SFRD in our re ferenc e model us¬ 
ing the SMHM relation from lGuo et a.1.1 (|2010f) . The solid 
red line is the SFRD that results from our approach of 
enforcing a monotonically increasing stellar mass in halo 
merger trees by retroactively suppressing star formation 
(see 42.21 a nd Appendix K|) . Obse rvations from the com¬ 
pilation of IBehroozi et all ( 201 3!) are also shown. The 
model SFRD lies above the data points at high redshift 
(z > 4) but below them at low redshifts {z < 2). 

Figure E^I b) shows the SFRD in our model using the 
redshi ft-depen dent SMHM function from IBehroozi et al.1 
(|2013l) . Since IBehroozi et~ahl (120131) forced their model 
to match the observed SFRD, one might expect that this 
SMHM function would also be consistent with the obser¬ 
vations. However, while our model SFRD (solid red line) 
is similar to the observations at high redshift (z > 6), it 
is still too low at lower redshifts. 

We have identified two significant effects that can raise 
the SFRDs in our models. First, additional star forma¬ 
tion is required to replace the stellar mass that is lost 
as massive stars reach the end of their lifetimes and die. 
We compute the additional star formation from t he sim ¬ 
ple ap proximation given in Eqn. (14) of Be hroozi et al.l 
(j2013r> for the fraction of mass lost from a single starburst 
as a function of age; this function is integrated over time 
to determine the ongoing star formation rate with recy¬ 
cled material. The resulting enhanced SFRD, shown in 
Figure [231 as the dashed blue line, is increased by factors 
ranging from 1.5 at z ~ 3 to 2.5 at z~0. This recycling 
effect is important and should be incorporated into fu¬ 
ture models, but it does not fully explain the difference 



















































Simulating Hubble Images With Galaxy Formation Models 


17 



o 

o 

CD 

-Q 

E 



Input mpieow 



0) 

o 


c 

o 

o 
£ 
0 ' 
"D 


20 


25 



30 20 

Input m F160W 


25 


30 



Fig. 19.— SExtractor number count (top row) and light (bottom row) detection efficiencies as a function of the input mpi6ow magnitudes 
(left) and Petrosian radii (right) of detected galaxies in the HUDF-depth simulated HST image, shown for each of the models explored in 
this paper. Note that the efficiency is lower for small galaxies because smaller galaxies are also fainter. 


between the computed SFRDs and the observations. 

The second potentially important effect is the ejection 
of stars from galaxies during merg ing. This is a crucial 
feature of the lBehroozi et all (1201311 calculation: the frac¬ 
tion of stars that escape from galaxies is an adjustable 
parameter. Those stars are then subsequently replaced 
by additional star formation. This parameter enables 
fitting both the observed SFRD and the galaxy masses 
in massive cluster-scale halos. Most of the intergalactic 
stars end up in clusters of galaxies and hence contribute 
to the diffuse intracluster light (ICL). 

To include the possibility of star formation enhance¬ 
ment through ICL, we adopt a more conservative proce¬ 
dure. Our standard approach with star formation sup¬ 


pression reduces the star formation rates in merger trees 
to ensure that the stellar mass of galaxies never de¬ 
creases. An alternative model is to allow all galaxies 
to form the stars required by the SMHM relation, and 
then to eject the excess mass (if any) to the ICL. The 
dot-dashed green lines in Figure [23] show the enhanced 
SFRDs that result from this pro cedure. With this ad - 
dition, the model SFRD for the iBehroozi et all ([2013T ) 
SMHM relation has been increased by another factor of 
2 and agrees reasonably well with the observations at all 
redshifts. 

Does this agreement imply that there must be a large 
fraction of stellar mass ejected to the ICL in order 
to match the observed SFRD? Probably not. In the 
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Fig. 20.— Detection efficiencies as in Fig. 1 191 but plotted as a func tion of redshift. We include only galaxies brighter than the median 
m Fi60W a t the stellar mass incompleteness limit, as shown in Fig. I18l 
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Fig. 21.— Apparent F160W mag nitude distribution from simulated HST images at HUDF depth for all models tested in this paper. 
Symbols are the same as in Fig. 1151 The biggest differences in the shapes of the simulated distributions come from the use of different 
stellar mass-halo mass relations (Models 1-4). Removing dust or metal content (Models 5 and 6) mostly shifts the distribution toward 
brighter magnitudes. 


IBehroozi et al.l (120131) computations, 30% of the stel¬ 
lar mass in the universe lies outside galaxies, with the 
ICL mass fraction being higher in clusters (P. Behroozi, 
private communication). Such a large fraction of ex- 
tragalactic stars may be in conflict with observations 
of the ICL, which su ggest ICL stellar mass fractions 
close r to 10% ie.g.. iMihos et ahl 120051 iZibetti et all 
20051: iK rick fc Bernsteinl 12007 : iGonzalez et al.l 120071 : 
Presotto et al.l 12014 iMontes fc Truiillol 120141 ). A model 
that relies on hiding much of the stellar mass outside 
of galaxies in order to enhance the SFRD may conflict 
as strongly with the observations as one that does not 


generate much intracluster light but has a lower SFRD. 

Instead of treating the inferred SFRD and ICL as 
known quantities, a better approach would be to apply 
the forward modeling methodology discussed in this pa¬ 
per: from simple theoretical models of galaxy masses 
and star formation rates, create simulated observations, 
including the SFRD and the ICL. Then the comparison 
between the models and the observations should be car¬ 
ried out entirely in the observational domain. It seems 
quite plausible that a consistent model can be created 
that fits the observed galaxy masses, observed SFRD in¬ 
dicators, and measured ICL properties; such a model will 
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Fig. 22. — Siz e functions from simulated HST images at HUDF depth for all models explored in this paper. Symbols are the same as 
those in Fig. eh Reasonably good agreement for larger objects is found with the real observations (black line) except for Model 7, which 
has larger galaxy sizes that do not scale with the halo size evolution. 
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Fig. 23.— Cosmic star formation rate density (SF RD) as a function of redshift. (a) Reference SMHM model frorii. (b) 
Redshift-dependent SMHM model from lBehroozi et alJ 1120131 1. The black dots show observed SFRDs compiled bv IBehroozi et all (l2013j l. 
The solid red line is the SFRD for our model with suppression of past star formation to create a monotonically increasing stellar mass in halo 
merger trees. The dashed blue line includes the additional star formation from instantaneous recycling of mass lost from shorter-lived more 
massive stars to maintain the stellar mass. The dot-dashed green line results from assuming that rather than suppressing star formation 
to enforce a monotonically increasing stellar mass, all halos form the stars required by the SMHM relation but the excess stellar mass is 
ejected from the halos to form intracluster light (ICL). The models including ICL fit the observations best but require that most stars 
reside in the ICL rather than in galaxies (see text for discussion). 

not necessarily demand a very large stellar mass in the 

ICL. 

5. DISCUSSION 
5.1. Summary of Results 

This paper demonstrates that a credible simulated uni¬ 
verse, consistent with deep HST images, can be con¬ 
structed using only a few basic assumptions and a small 
number of parameters. The simplicity of our model, in 


contrast with the complexity of most other models of 
galaxy evolution, enables robust comparisons with the 
observed universe. A key element of this comparison is 
to project our models of evolving galaxy populations into 
the observational domain and to extract measurements 
from the simulations using the same tools as used to an¬ 
alyze the real images (“forward modeling”). Since our 
simulated HST images include all relevant cosmological, 
instrumental, and selection effects, they can be compared 
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directly with real HST images. In order to simulate 
the dumpiness and internal structure of galaxies, we cut 
out images of SDSS galaxies as templates for our model 
galaxies (instead of the commonly used smooth Sersic 
profiles), with their fluxes and sizes rescaled to reflect 
the properties of model galaxies in our simulations. 

The key assumption of our theoretical modeling is that 
most of the information needed for a first order match be¬ 
tween the simulated and observed universes is already en¬ 
coded in the gravitational dynamics of dark matter halos, 
as derived from cosmological IV-body simulations. In this 
semi-empirical approach, we express the stellar mass and 
size of a model galaxy as a function of the mass and size 
of its host halo. This effectively avoids any detailed mod¬ 
eling of complex baryonic physical processes, in contrast 
to semi-analytical and hydrodynamical simulations, with 
their larger number of assumptions and parameters. Our 
semi-empirical models are based on statistical matching 
between the distribution of halo mass and size measured 
in the Wbody simulation to the distribution of stellar 
mass and size for real galaxies at z = 0, thus guarantee¬ 
ing a good match between the present-day simulated and 
observed galaxy populations. For the evolution of galaxy 
populations, we adopt several simple models for the red- 
shift dependence of the stellar mass-halo mass (SMHM) 
relation. We also show that the evolution of the lumi¬ 
nosities and spectra of galaxies can be modeled, to good 
approximation, using the star formation histories implied 
by the growth of the stellar mass of each galaxy along the 
merger tree of its host halo. The galaxy spectra are de¬ 
termined by convolving the star formation histories with 
stellar population synthesis models. 

The analysis of simulated HST images is a powerful 
tool for comparing the predictions from different mod¬ 
els of galaxy evolution. By comparing directly in the 
observational domain, we can decide which choices of 
parameters make an observable difference and which do 
not. Our reference model assumes a non-evolving SMHM 
relation and non-evolving solar metallicity and dust con¬ 
tent. Despite these simple assumptions, it provides a 
good match to real HST images, particularly when com¬ 
paring the luminosity and size distributions. Including 
plausible redshift-dependence in the SMHM relation does 
not radically alter the simulated HST images. In con¬ 
trast, a much less plausible, linear SMHM relation, with 
the stellar mass proportional to the halo mass, leads to a 
far greater abundance of compact, high-luminosity galax¬ 
ies in low-mass halos at high redshift. From this and our 
other models, we conclude that the stellar mass efficiency 
per unit halo mass should indeed peak at stellar masses 
around M* and decrease for smaller and higher masses. 
Metallicity and, especially, dust have a strong effect on 
the observed galaxy properties, with their luminosities 
increasing by up to 0.7 mag when removing dust or met¬ 
als from the galaxies. From the analysis of galaxy sizes, 
we conclude that galaxies must be smaller in the past. 
A simple linear scaling between galaxy and halo sizes 
provides a good match between the size distributions of 
simulated and real HST images. If there were no evolu¬ 
tion in the sizes of galaxies, the deep HST images would 
appear almost empty. 

We also find that the measured values of size and lu¬ 
minosity of galaxies in the simulated images are strongly 
biased. SExtractor underestimates the luminosities and 


sizes of simulated galaxies especially around the magni¬ 
tude detection limit, as the extended low surface bright¬ 
ness components of galaxies are lost in the image back¬ 
ground noise. These biases vary among the models, 
depending primarily on the size-luminosity relation of 
galaxies. For example, our model with a linear SMHM re¬ 
lation predicts more light for small, faint galaxies, which 
makes them easier to detect compared with our reference 
model. 

The detection efficiency of galaxies, measured by com¬ 
paring the number and luminosities of detected galaxies 
in the simulated images to those for all model galaxies, 
also reveals some interesting results. For the reference 
model, the number detection efficiency declines slowly 
with magnitude as objects get fainter to ~ 80% just 
above the detection limit and it then drops more rapidly. 
This behavior is similar at low and high redshifts, since 
it depends only on whether galaxies are bright enough to 
pass the magnitude cut. The fraction of light recovered 
from galaxies, on the other hand, falls considerably faster 
with magnitude; it is ^60-70% right above the detection 
limit and then drops sharply. The reason for this is that 
SExtractor not only fails to detect faint galaxies but 
also underestimates the fluxes of the ones it does detect. 
Note that although the model with a linear SMHM rela¬ 
tion behaves similarly as a function of magnitude; it has 
a much higher detection efficiency at small sizes, since it 
assigns more luminosity to the small, high-redshift galax¬ 
ies, making them easier to detect. The fraction of missing 
light may be underestimated in our simulations because 
they truncate the luminosity profiles of galaxies at two 
Petrosian radii. 

The redshift dependence of the detection efficiencies is 
also interesting but requires some additional considera¬ 
tions. The numerator of the efficiency (i.e., the number 
of detected galaxies or the total light detected) is mea¬ 
sured with reasonable accuracy. However, the denomina¬ 
tor (the total number or luminosity of galaxies) is not well 
determined because our images are missing galaxies with 
masses below the resolution of the milli-Millennium simu¬ 
lation. We address this issue by fitting the power-law tail 
of the luminosity function, and estimate the missing light 
by extrapolating the fitting function toward zero lumi¬ 
nosity. Our measurements for the reference model show 
that the fitted power-law slope becomes steeper at higher 
redshift, changing from a ~ —1 at z < 1 to a = —1.75 
at 2 = 7 — 8. Since the slopes are shallower than the di¬ 
vergence value a = —2, the extrapolated missing light is 
always finite and the fraction of missing light is modest 
in our case. Indeed, the light detection efficiency drops 
from ~90% at small redshifts to only ~50% at 2 = 7 — 8. 
These values depend sensitively on the slope of the lumi¬ 
nosity function, which in turn depends on the dark mat¬ 
ter halo mass distribution. Since the power-law tail of 
this mass distribution has a slope of a = —2 (confirmed 
both by theory and simulations), it would not be diffi¬ 
cult to find models with suitable combinations of SMHM 
relations and mass-to-light ratios that predict a power- 
law slope of the luminosity function closer to a = —2 or 
even steeper, resulting in very small detection efficien¬ 
cies. We demonstrated this by our model with a linear 
SMHM relation, which presents a lower light efficiency 
than models with more plausible SMHM relations. The 
total number of galaxies, on the other hand, diverges for 
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slopes steeper than the limit a < — 1. Since our simu¬ 
lated universe presents steeper slopes, we conclude that 
the number counts detection efficiency as a function of 
redshift cannot be accurately measured using the milli- 
Millennium simulation at most redshift. 

5.2. Future Directions 

The main goal of this paper was to develop some first- 
generation simulations of deep HST images using semi- 
empirical models of evolving galaxy populations. In this 
spirit, we have tried to strike a balance between mod¬ 
els that are realistic enough for meaningful conclusions 
and models that are simple enough for computational ef¬ 
ficiency and ease of interpretation. Indeed, we regard the 
simplicity of our models, with relatively few assumptions 
and parameters, as a definite virtue in a preliminary ex¬ 
ploration such as this. Now that this initial analysis has 
been successfully concluded, it is appropriate to consider 
how our approach can be enhanced for greater physical 
realism and accuracy. In the remainder of this section, we 
outline several directions for future studies of this type. 

A relatively straightforward enhancement of our mod¬ 
els would be to base them on dark-nratter simula¬ 
tions with larger comoving volumes and smalle r particle 
mass es, such as th e Millennium II dBovlan-Kolchin et all 
120091) or Bolshoi (IKvplin et all 1201 ill simulations. The 
larger volume would permit a more accurate simulation 
of both the small and large-scale distribution of galaxies 
in deep images than is possible with the milli-Millennium 
simulation. In particular, it should be possible to lay 
down light cones through the simulation volume to com¬ 
pute galaxy distributions rather than relying on sampling 
the halo distribution as we do in this paper. Such simu¬ 
lations, in addition to providing direct estimates of cor¬ 
relation functions and other clustering statistics on large 
scales, will provide estimates of the variance in counts 
and other properties of galaxies in smaller fields. The 
higher mass resolution of larger dark matter simulations 
will enable the modeling of galaxy populations to lower 
masses and luminosities. This is especially important 
for the simulation of JWST images, which will reach 
much fainter magnitudes than the HST images consid¬ 
ered here. The extension of model galaxy populations 
to lower masses will also improve estimates of detection 
efficiencies for the counts and light of galaxies, which de¬ 
pend on extrapolations of luminosity functions below the 
detection thresholds. 

It would be also worthwhile to make simulations with 
a larger suite of SMHM relations. The resulting sim¬ 
ulated HST images could then be compared with ob¬ 
served images and a goodness of fit assigned to each of 


these SMHM relations (by e.g. maximum likelihood). 
The requirement that the models satisfy constraints on 
the SFRD and ICL (§4) could be added to the pro¬ 
cedure at this stage. In this way, the “best” SMHM 
relation and confidence regions around it could be de¬ 
termined. As a byproduct, such statistical tests would 
indicate how much useful information HST images re¬ 
ally contain about galaxy formation and evolution. In 
particular, we would like to know whether the evolu¬ 
tion of the SMHM relation can be pinned down uniquely 
or whether substantially different SMHM relations lead 
to similarly good fits in comparison with HST observa¬ 
tions. This issue has important theoretical implications; 
different SMHM relations presumably reflect differences 
in the physical processes in galaxy formation and evo¬ 
lution. This issue also has practical implications, be¬ 
cause simulations like those presented here could be used 
to guide future observing strategies, with different tele¬ 
scopes, cameras, and filters, to determine the SMHM 
relation most efficiently and robustly and to resolve am¬ 
biguities in the models. 

Another improvement would be to allow the metal and 
dust content of model galaxies to evolve with redshift. 
Our results, based on non-evolving metal and dust con¬ 
tent, show that differences in these quantities are read¬ 
ily apparent in simulated HST images. A physically 
plausible scheme for evolving the metal and dust con¬ 
tent should be based on the star formation and thus the 
metal enrichment history of each galaxy in the model. 
Unfortunately, the stellar metallicity and ISM dust opti¬ 
cal depth also depend on the evolution of the ISM of each 
galaxy, including inflows and outflows, and these can¬ 
not be specified without introducing additional assump¬ 
tions and parameters. This is why we have neglected the 
evolution of the metal and dust content in the present 
exploratory study. Nevertheless, such effects could and 
probably should be included in future studies, provided 
one is willing to pay the price of extra assumptions and 
complexity. 
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APPENDIX 

A. ENFORCING A SELF-CONSISTENT, MONOTONICALLY INCREASING STELLAR MASS IN THE MERGER TREE 

We have developed a self-consistent and recursive method that ensures a monotonic increase in the stellar mass 
content of dark matter halos across their merger trees. This retroactive suppression mechanism functions in the 
following way: when a halo is measured to decrease its current M s value, we suppress accordingly the value of M s 
in its progenitor halos located in the closest previous time steps, which consequently causes scatter away from the 
one-to-one SMHM relation. First, we start with a root halo picked from the collection {hi}lZi of all N halos found 
in the last (fcth) time step in the simulation (most likely to be at z = 0), and compute M s = M s (Mh a io) for all halos 
in its merger tree using the one-to-one SMHM relation. Then we pick its progenitors located in the previous (k — l)th 
time step and verify that EM SjProg < M s>root . If not, we reduce the stellar masses of these progenitors (proportionally 
to their masses), so that EM S)Prog = M Sj desc holds. We continue decreasing the M s values in halos of all the closest 
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previous time steps until reaching the jth time step (j < k ), where EM SiProg < M SiIoot holds for all progenitor halos 
located at the (j — l)th time step. This process is repeated recursively and backward in time throughout the merger 
tree of the root halo, with each one of its progenitors playing the role of the root halo, that is, each progenitor goes 
later under the same process that modifies M s values along its own (smaller) merger tree. All the previous steps are 
then repeated for the next root halo in {/i,;}*E) v until the whole dark matter simulation is traversed. Evidently, the 
optimal case would be to have the fcth time step located (if available) at sometime far enough in the future, so that 
all values of M s at z = 0 are modified and stabilized. 

B. GALAXY IMAGE CUTOUTS FROM SDSS 

This appendix describes the procedure used to select SDSS galaxy image cutouts that represent model galaxies on 
our simulated HST images. __ __ 

As a repository of galaxy images, we use the SDSS Data Release 10 (lYork et alj|2000t lAhn et aLll2014H . This data 
archive (a MS-SQL Server database available online via CasJobtQ) is suitable for our needs, as it provides imaging and 
spectroscopy of ~ 10 4 5 6 7 galaxies as well as many other d erived physical pro perties. 

We use in particular the Main Galaxy Sample (MGS. lStrauss et aj]l2002 ). which spans a rich variety of galaxy types 
over a sky area of 7930 deg 2 . These galaxies form a flux-limited sample, with an r-band Petrosian apparent magnitude 
limit of m r < 17.77. We define our sample in the redshift range of [ 21 , 22 ] = [0.01,0.2], with the peak of the redshift 
distribution located at 2 ~ 0.1. We consider clean science primary galaxies only on calibrated images having the 
photometric status flag. We discard suspicious detections and objects with deblending and interpolation problems, as 
well as objects whose spectral line measurements and properties are labeled as unreliable^- Our main sample contains 
~400,000 galaxies. 

The properties available for the MGS that we use in this paper are the photometry in the ugriz bands, in particular 
model magnitudes for computing the colors (iStouehton et ali [2002h . A measure of the galaxy size is given by the 
(redshift insensitive) r-band Petrosian radius Rv> (jPetrosianlll976l l and its less noisy proxy R 50 (the radius encircling 
50% of the Petrosian flux). The galaxy stellar mass M s is also available; it was obtained from spectral analysis at 
MPA and JH10- as detailed bv iKauffmann et all (12003H . 

We extract the galaxy cutouts from 50 TB of 5-band, 2048x1489 pixel SDSS image frames us ing par allel processing- 
in the Data-ScopeQ a computing facility designed specifically for data-int ensive ap plicatio ns dSzalav et al.ll2012il . As a 
first step, we denoise the images using the anisotropic diffusion method (jPerona fc Maliklfl990i l, which is particularly 
useful for the more noisy u and 2 bands. This leaves the galaxy outline and internal structure intact, since the 
smoothing is done in the direction locally parallel to the edges or borders. Multi-band images are registered to the 
same reference frame using cubic spline interpolation. 

We extract individual galaxy cutouts from the original images, considering only what is included inside an aperture 
of radius 2Rp centered at each galaxy. To avoid the influence of neighboring objects, we do not use galaxies where the 
sum of fluxes from all neighbors that overlap the aperture exceeds 1% of the galaxy flux. This leaves us with a total 
of ~270,000 galaxies in the repository. 

We select from the SDSS image database the galaxy that is the closest neighbor to the model galaxy in the multi¬ 
dimensional space of galaxy properties. For this paper, we use two sets of parameters to select matching SDSS 
galaxies: 

1. u — r color and log M s 

2. u — r color, log M s and logi? 5 o 

Other schemes using additional parameters were tried with very similar results. Each property is mean-subtracted 
and normalized by the standard deviation before doing the matching. The u — r color tracks the flux difference before 
and after the 4000A break, and is a good indicator of age. To compare photometry, the (evolved) rest-frame spectrum 
of the model galaxy is redshifted to the median of the SDSS redshift distribution (z = 0.1), and then convolved with 
the u and r filters. Using the properties from a simulated universe based on our reference model ([j3]), we find that 
the median of the rms error over all properties derived from the second matching scheme (logrms = —0.85) is bigger 
than that of the first scheme (logrms = —1.18), since there is one more variable to match (log-Rso). However, the 
individual median rms error for logi?so in the second matching scheme is ~ 1 dex smaller than that of the first. A good 
compromise is to use the first matching scheme to get a smaller overall rms error, but then to manually rescale the 
SDSS galaxy size to match the model galaxy size as described in 112.61 This approach is used as part of our reference 
model as well. 
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